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Abstract 

V-^ \ This work presents new methods and algorithms for tracking the shape 

. and trajectory of moving reflecting obstacles with broken rays, or rays 

00 ' reflecting at an obstacle. While in tomography the focus of the recon- 

, struction method is to recover the velocity structure of the domain, the 

shape and trajectory reconstruction procedure directly finds the shape 
m^j ' and trajectory of the obstacle. The physical signal carrier for this innova- 

tive method are ultrasonic beams. When the speed of sound is constant, 
\^ • the rays are straight line segments and the shape and trajectory of moving 

C/3 ' objects will be reconstructed with methods based on the travel time equa- 

O . tion and ellipsoid geometry. For variable speed of sound, we start with the 

eikonal equation and a system of differential equations that has its origins 
in acoustics and seismology. In this case, the rays are curves in R 3 and we 
develop algorithms for shape and trajectory tracking based on the numer- 
ical solution of these equations. We present methods and algorithms for 
£N| , shape and trajectory tracking of moving obstacles with reflected rays when 

C**) ■ the location of the receiver of the reflected ray is not known in advance. 

' The shape and trajectory tracking method is very efficient because it is 

not necessary for the reflected signal to traverse the whole domain or the 
same path back to the transmitter. It could be received close to the point 
of reflection or far away from the transmitter. This optimizes the energy 
spent by transmitters for tracking the object, reduces signal attenuation 
and improves image resolution. These are some of the key advantages for 
using reflected rays in designing this new class of imaging and tracking 
■ methods and algorithms. Another key advantage and property of shape 

' and trajectory tracking is that it works for moving obstacles and thrives 

on the movement and dynamics of obstacle, transmitters and receivers. 
It is a safe and secure method. We also present algorithms for tracking 
the shape and trajectory of absorbing obstacles. The new methods and 
algorithms for shape and trajectory tracking enable new applications and 
an application to one-hop Internet routing is presented. 



1 Introduction 

Consider a bounded convex domain ilo £ R 3 with one moving obstacle fii(t) C 
ilo, where t £ R. Let f2(t) = firj\^i(i) where Qi(t) is a convex set in R 3 with 
a smooth boundary. Let f2i(t) C fio- Suppose that for all t we are given all 
integrals J f(l)dl = C 1 where 7 are broken rays in f2(i) such that for each 
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Figure 1: Ellipsoid Geometry 



point P G dfli(t) there is at least one broken ray 7 reflecting at P. A broken 
ray is a ray reflecting at the obstacle and starting and ending at the observation 
boundary dflo. Let f(x) — v ^ x ^ > in f2. Then as we know C 7 correspond 
to signal travel times in a medium with speed of sound «n (x) . The shape and 
trajectory reconstruction problem is to find d£li(t) given the set C 1 (t) or the 
functions Ct("f) where 7 e Ct(t). This problem has two separate cases. The 
first case corresponds to finite travel times of signals reflected from a reflecting 
obstacle. The second case is for infinite travel times that correspond to signals 
that are not received at any receiver. One example of such signals can be 
observed in the case of obstacles that absorb signals instead of reflecting the 
signals. In this second case, I will give an algorithm for reconstructing the 
shape and trajectory of the objects as the intersection of beams absorbed by 
the obstacles. In the first case of finite travel times and constant speed of sound, 
the shape and trajectory of moving objects will be reconstructed with methods 
based on the travel time equation and ellipsoid geometry. I then solve the shape 
and trajectory reconstruction problem for finite travel times and variable speed 
of sound with methods derived from the eikonal equation. 



2 Ellipsoid geometry in M 2 

Let E be an ellipse in M 2 with focci F\ and F2. For any point P <S E, as 
illustrated in Fig. 1, the defining geometric property of the ellipse gives 

+ \F 2 P\ = 2a. (I) 

Choose a Cartesian coordinate system such that F\ = (— c, 0), F2 = (c, 0) 
and let E have major axis with length 2a, where a > c > 0. Then E is described 
by the equation 

~2+tl = ^ ( 2 ) 
a 1 b l 

where b 2 = a 2 — c 2 and 2b is the length of E'a minor axis. 

Let |i*iP| = r and P = (x,y). From the analytical geometry of conic 
sections, 

r = a H — x. (3) 
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Indeed, observe that because larl < a and < - < 1, it follows that a + -x > 0. 
Therefore, 

r 2 = |FiP| 2 = x 2 + y 2 + 2xc + c 2 = 

b 2 x 2 + b 2 v 2 9 

= ryr- — + 2xc + c 2 = 

b l 

b 2 x 2 + (a 2 - c 2 )y 2 „ , 
b l 

b 2 x 2 a 2 y 2 2 y 2 

= — + ^ + 2 - + c (!-^) = 

b 2 x 2 + a 2 y 2 



b 2 

a 2 b 2 



2xc + c^ — 

2 

2xc + c 2 %- 



b 2 

„2 



= a 2 + 2xc + c 2 ^r 



(a+-x) 2 



implies that r = a + -x 

If F\P makes an angle 9, < 9 < 2t:, with the positive x-axis of the 
coordinate system then 



x = r cos 9 — c 
y = r sin 9 



(4) 



Combining 3 and 4 gives 



or 



c 

x = r cos 9 — c = (a H — x) cos 9 — c 
a 



a cos a — c 

x 



(5) 



Therefore, 



1 - ^ cos 9 

a 

c , a cos 6 — c . 

r = a + -{- -) 6 

a 1 — r cos ^ 



and 



a cos 9 — c 
X = 1 - ^cos6» 

a 

. c a cos 9 — c . . . 

y = (o+ -(- E 7 sinfi 

a 1 — - cos 9 



(7) 



3 Ellipsoid geometry in R 3 

Let E be an ellipsoid of revolution in R 3 with focci F\ and F2 . For any point 
P € E 1 , as illustrated in Fig. 1, the defining geometric property of the ellipsoid 
gives 

\FiP\ + \F 2 P\ = 2a. (8) 
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Choose a Cartesian coordinate system such that F\ = (— c, 0, 0), F2 = 
(c, 0,0) and let E have major axis with length 2a, where a > c > 0. Then 
E is described by the equation 



9 2 9 

x y z 

where b 2 = a 2 — c 2 and 2b is the length of E's minor axis. 

Let \F\P\ = r and P = (x,y,z). From the analytical geometry of conic 
sections, 

r = a+ -x. (10) 
a 

If F\P makes an angle tf>, < (j> < tt, with the positive z-axis, and the 
projection of F\~P to the x-y plane makes an angle 6, < 6 < 2ir, with the 
positive 2;-axis of the coordinate system then 

x = r sin <f> cos 6 — c 

y = r sin <p sin 9 (11) 
z = r cos 4>- 

Combining 3 and 11 gives 

x = r sin <j> cos 9 — c = (a H — x) sin <fi cos 9 — c 

a 



ox 



a sin m cos ( 



■ sin d> cos 9 



Therefore, 



and 



1 — r sin (A cos ( 



(12) 



c . a sin <p cost* — c . . 
a+ a ( l-^sin ( icos0 ) (13) 



a sm cos ( 



. c. asm cost/ — c.. . , . . , , 

y=(a+-{- ^))sm0sm6» (14) 

a 1 — - sm <b cos 9 y 

a t 

, c , a sin d> cos 9 — c . . , 

z = a + - . cos0. 

a 1 — - sm 6 cos 6* 



4 Reconstruction for constant speed of sound 

Let L G 9r2o and P G 90i and send an ultrasonic beam from L to P. Let 
be reflected at P and denote the reflected beam as P$ , where S G <9£!o is 
the intersection of the reflected beam and dflo. In addition, place synchronized 
clocks at L and S and measure the time t that it takes the ultrasonic beam to 
travel the distance \LP\ + \PS\. The time of flight t can be measured by other 
methods as discussed in section 10. Denote the speed of the beam in f2 as vq. 
Then 

\LP\ + \PS\ = tv Q = 2a (15) 



4 



where a > 0. When L and S are fixed, equation 15 for the position of point 
P is the equation of an ellipse in R 2 or the equation of an ellipsoid of revolution 
in M 3 . I will refer to equation 15 as the travel time equation. 

Let \LS\ = 2c, where c > 0. The triangle inequality implies a > c. Again, 
consider a Cartesian coordinate system with center O such that L = (— c, 0,0) 
and S = (c, 0, 0). Measure the angle <f>, < <j> < tt, that LP makes with 
the positive z-axis of a coordinate system centered at L and parallel to the 
coordinate system with center at O, and measure the angle 0, < < 2n, that 
the projection of L~P to the x-y plane makes with the positive x-axis of the 
coordinate system centered at L. Then the coordinates of P with respect to this 
coordinate system centered at O are given by 14. 

A transmitter T at point L C dilo, emits signals Sk, < k, that are received 
by a receiver R at point S C 90o- The signals may be reflected or absorbed 
by fii on their way from T to R. The time tk that each signal Sk travels from 
T to R can be measured and is considered a known quantity. The shape and 
trajectory reconstruction problem at hand is how to reconstruct as a function 
of time the boundary surfaces dQ\ from the set {tk}- This problem can be 
reduced to the problem of reconstructing the shape of the obstacle dfti from 
measurements {tk} taken during a fixed sampling period r during which the 
position of the obstacle will not change significantly. The trajectory of the 
obstacle is reconstructed by successively reconstructing the shape of the obstacle 
9J7i(ri) for successive discrete sampling periods t,, where 1 < i < oo. From 
now on, the measurement set {tk} will refer to a specific sampling period r. 

This section considers the shape and trajectory reconstruction problem when 
each signal Sk is a beam, such as an ultrasonic beam, that travels along short 
linear paths and is reflected from a single point Pk in dQ,\. In this restricted 
version of the problem, we consider the case when the path of each signal Sk is 
piece- wise linear with a single reflection point Pk and when each Sk travels with 
constant speed along the two sides TkPk and PkRk of ATkPkRk in K 3 - Then 
there is an elegant geometric solution to our reconstruction problem: 

The travel time equation allows the reconstruction with finite resolution of 
dfli by sending from T = L multiple beams and measuring the arrival times of 
these beams at R = S. When sending beams {bk} at known angles <fik and 0k 
and measuring the corresponding travel time tk of the beams bk we can compute 
a k = ^f 51 and then Xk, yk and Zk by formulas 14: 



a k sin <pk cos t) k 
Xk = V" 



- sin 4> k cos 9 k 

u -k 

c a k sin<f>kCos6k - c.. . . 

Vk = (ak H h -. — : — 7 — ))sin0 fe sin0 fc (16) 

a k 1- ^Sin0 fc cos0fc ' 

c a k sm(j>kCOs9k -c.. 

Zk = (Ofc H (- c . — ))cos0 fc . 

Ofe 1 - — sin <pk cos k 

The angles 4>k and 0k can be generated according to different methods for 
scanning Q,\. For example, let N be a large positive integer that specifies the 
level of resolution for <pk and 0k- Then we can set 

fa = jjl,0 < K N 
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and 

6 k = —m, < m < N 

Alternatively, Ok and <f>k can be generated by other methods. For example, 
9k and <pk can be generated randomly. The transmitter can send rays in different 
directions in parallel at the same time or in a sequence one ray at a time. 

The method described in this section implies that the obstacle's boundary 
acts as a mirror. One can then ask why it is necessary to use reflection and el- 
lipsoid geometery in order to reconstruct the shape and trajectory of the mirror. 
The reflected signal may be received at the transmitter and the reflection point 
may be reconstructed by computing the distance along the transmitted ray for 
half of the time of flight. When the signal is reflected so that that angle of inci- 
dence equals the angle of reflection and the shape of the obstacle is irregular, the 
reflected signal will not necessarily go back to the transmitter. Moreover, the 
reflecting obstacle oscillates or vibrates due to mechanical or other reasons. The 
most basic mathematical model that we can choose to model such oscillations 
is the simple harmonic oscillator. 

d 2 9 

-jjp- + awnBi = 

where a is a constant coefficient and Qi is the angle between the incoming 
ray and the surface normal at reflection point P. When the tangent plane to P 
does not oscillate, 9i does not oscillate and the reflected ray may be received 
close to the transmitter. However, when the tangent plane to P oscillates, the 
reflected ray is received at varying distances from the transmitter. For a system 
of coupled oscillators 

9 1 = f 1 (t,e 1 ,e 2 ,...,e n ) (17) 

: (18) 

On = fn(t, 01, 9 2 , On) (19) 

the angles of reflection can be chaotic. The movemements of the tangent 
plane at the reflection point P can be due to turbulence. Therefore, it is hard 
to control the angle of incidence and the angle of reflection. 

Reconstruction via ellipsoid geometry allows reconstruction of P with a 
transmitter and a receiver. In practice, it is not clear where to position the 
receiver. One approach to solving this problem is to create a distributed sys- 
tem of receivers. When the receivers are located sufficiently close to each other 
on the observation boundary, every reflected signal will be received by some 
receiver. This requires the placement of a lot of receivers on the observation 
boundary. One approach to minimizing the number of receivers is to consider 
a system of receivers moving on the boundary. For example, in a flow on the 
boundary with speed Vf a receiver moving with the flow will cover the boundary 
in time 

i(dn ) 

Vf 

where Z(<9f2o) is the length of dflo. N receivers that are distributed uniformly 
and move in the same direction can cover the boundary in time 
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When t c < i m inj where t m - ln is the minimum time for which the transmitter is 
continously sending a given unit of information from the signal, then all reflected 
signals will be detected. This method for minimizing the number of receivers 
on the observation boundary can be considered a distributed algorithm: the 
receivers working as a team move on the observation boundary in order to 
minimize the number of receivers. Another way of looking at this approach 
is to consider it as a method of converting the energy of the boundary into a 
function that minimizes the number of receivers. Let M = Nm be the total 
mass of the stationary receivers needed for detecting one bit of information. Let 
the kinetic energy of a receiver on the boundary be E/.- Then the velocity of 
one receiver is 

2E k l2E k N 



M 

Substituting this expression into 20 gives 



(21) 



For t c = i m in we need at least 

N= iim^ (22) 

'nun 

receivers uniformly distributed along the observation boundary and moving 
with speed v m in the same direction, clockwise or anti-clockwise, along the 
boundary. 

In this work we consider a distributed system of transmitters and receivers 
located at the observation boundary such that each transmitted signal is received 
by one receiver. 



5 Reconstruction for variable speed of sound 

Consider a wave, or signal, described by the wave equation 

u tt - c 2 (x,y, z)Au = (23) 
where c(x, y, z) > is the variable speed of sound in Q(t) = f2o\Oi(i) and 

«|afii(t) = 



for t > and Qi(t) C Oo C R , where fii(£) is a moving convex obstacle 
with a smooth boundary. 

We consider an environment without caustics and look for solutions of the 
form 

00 e iu(W(x,y,z)-t) 

u(x,y,z,t) = ^2Aj(x,y,z) — (24) 

3=0 ^ 

where the eikonal function 

W(x, y, z) = const 
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defines a surface of constant phase. These solutions of the ray equation are 
called rays or ray solutions. 

Plugging 24 into the wave equation and in the resulting series equating 
the coefficients of the powers of cj on both sides of the equation gives for the 
coefficients of the highest power uj 2 

u 2 A (x, y, z) = uj 2 A (x, y, z)c 2 {x, y, z)(VW(x, y, z)f (25) 

Setting to the coefficient of oj 2 Aq(x, y, z) in equation 25 leads to the eikonal 
equation 

(VW(x,y,z)) 2 = 1 (26) 
c(x,y,zy 

The eikonal equation will serve as the starting point for deriving equations 
for tracing the ray path. 

We follow [E, JG, R] to present the derivation of a system of first order 
ordinary differential equations from [JG] for tracing the path of the ray. The 
initial values for solving this system are the ray's initial position and take-off 
angles. We then apply the numerical solution of this system of ray tracing 
equations for the development of a novel algorithm for shape and trajectory 
tracking of obstacles. 

Denote by s the length of the ray from its initial position to point (x, y, z) 
along the ray path. Then 

^ = c(x(s), y(s),z(s))VW(x(s), y(s),z(s)) (27) 
as 

because by the eikonal equation c(x(s),y(s), z(s))VW(x(s),y(s), z(s)) is a 
unit vector normal to the wave front and tangent to the ray path. i?(s) = 
(x(s), y(s), z(s)) is the position vector of the ray path. 

From equation 24 it follows that W(x(s), y(s), z(s)) has units of time and rep- 
resents the traveltime of the wave from its origin to point i?(s) = (x(s), y(s), z(s)). 
Then because speed is the derivative of distance with respect to time 

^'•'■' W) = c(l( ,J, )t , ( , )r w.). rt .).,(. )) w 

The function U(x(s), y(s), z(s)) = ^jy^j^r^ is reciprocal of the speed 
and is referred to as the slowness of the wave. 

Taking the gradient of both sides of 28 and swapping the operators d and V 
gives equation 

dv^(,( s ) ; ,( g ) ; K g )) =w(a;(s)i ^ )iZ(s)) (29) 

Combining 29 and 27 leads to 

^(U(x(s), y(s), z(*))^) = Vtf (x(*), y(s), z(s)) (30) 
The unit vector m t ne direction of wave propagation is 



ds 

dR(s) 



(sin cj)(s) cos 0(s), sin </>(s) sin 8(s), cos 4>{s)) (31) 
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where R(s) is the position vector of the ray and where </> is the incidence 
angle between the ray and the z axis and 9 is the azimuth angle between the 
projection of the ray vector to the xy plane and the positive x axis. 

Combining equation 31 and 30 leads to the following equations 

= U cos 4> cos Q-7- — U sin (f> sin 6^f- + sin cf> cos 9-^— (32) 
dx as as as 

= U cos 6 sin 9^- + U sin 6 cos 9^- + sin 4> sin 6^- (33) 
dy as as ds 

dU rr ■ J<t> ± 9U 

— — = -(7sm0— + cos0— — (34) 
dz as as 

This system of three equations can be simplified to a system of two equations 
that does not include 

drh r-naW) QU . n dU . sin 6 dU . „. 



ds U dx dy U dz 

d6 1 , n dU . n dU, 

— = — — cos0— -sin0— 36 
as U sm<p ay dx 

In addition, we know that 

dx(s^ 



ds 
dy(s) 
ds 



- sin cf)(s) cos 9(s) (37) 

= sin cj>(s) sin 0(s) (38) 

^=cos0(s) (39) 
as 

Our data includes time of flight and we will transform the above five equa- 
tions to have time instead of arc length derivatives. This is accomplished by 
applying the following equations: 

d9 ^d9 d<f> ydrf) dx jj^ x dy jjdy dz ^dz 
ds dt ' ds dt 1 ds dt ' ds dt 1 ds dt 
dU__ _J_dcdU_ _J_dcdU_ 1 dc 

dx c(x) 2 dx' dy c(x) 2 dy' dz c{x) 2 dz 

As a result, we arrive at the equations for the Shooting Method for two-point 
seismic ray tracing [JG, R, Y]. 

dx 

— = c(x, y, z) sin cos (42) 
dy 

— = c(x, y, z) sin sin (43) 

dz 

— = c(x, y, z) cos <j> (44) 
dd> dc dc dc 

— = — cosd>(— cos6 + — sin#) + — sinei (45) 
dt dx dy dz 

d9 1 ,dc . . dc n . ,,„. 

-(— sinfl- — cosfl) (46) 



dt sin^ dx dy 
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This system of equations has wide applications in seismology and is used in 
algorithms for seismic ray tracing [SK]. 

The system of equations 42 is for a Cartesian coordinate system where R(t) = 
(x(t),y(t), z(t)) is the ray position vector, <p(t) is the incident angle of the ray 
direction vector with the z axis and 6(t) is the azimuth angle that the projection 
of the ray direction vector makes with the positive x axis . 

We consider that the speed of sound c(x, y, z) is known inside f2. We will 
reconstruct the position of the reflection point P given the positions of the 
transmitter and receiver and the zenith and azimuth angles of the transmitted 
ray with respect to a Cartesian coordinate system centered at the transmitter. 

We know the angle and initial speed of wave propagation for the transmitted 
signal at time t = i.e. we know the velocity of the ray at t = 0. We can also 
measure and know the initial position (x(0), 2/(0), z(0)). Then, knowing the 
initial conditions, we can find numerically the signal path from the transmitter 
L through the reflection point P for a given travel time Tk- Let the travel time 
from the transmitter L to the reflection point P be Tk and the total travel time 
from the transmitter L to the receiver S through the reflection poin P be tk 
where tk > Tk- 

From symmetry, we can imagine that the signal received at the receiver S 
is transmitted from the receiver S and that it arrives at the reflection point P 
in time t = tk — Tk because its travel time is tk — Tk- The ray path from S to 
P is again described by the above system of equations for seismic ray tracing, 
however we do not know the initial angles for the signal from the receiver S. 

In order to reconstruct the intersection point Pk of rays starting from the 
transmitter and receiver for a given data point 

B k = (xi,yi,zi,x r ,y r ,z r ,(f>k,9k,tk,£k) 

containing the measured values for the signal, we step through a discrete set of 
values for Tk , where <f>k and 9k are the zenith and azimuth angles of transmission, 
xi, yi, zi are the coordinates of the transmitter endpoint of the ray, and x r , y r , 
z r are the coordinates of the receiver endpoint of the ray, tk is the time of flight 
for the signal and is a frequency of the signal. Starting to trace a ray path 
from L, at each time step, each new point on the path is a candidate reflection 
point P. At this stage, this is basic initial value ray tracing. Then we apply 
the two point seismic ray tracing shooting method to see if we can reach from 
receiver S candidate point P for the remainder of our time budget. If we can 
then we have found the reflection point P for this data point. Otherwise, we 
continue with the next time step of the initial value ray tracing and check the 
next candidate point P. We repeat this until we find a reflection point P or 
exhaust our time budget tk- 
For example, 

_ t k n s 
Tks ~ N r 

where N r is an integer that specifies the time-step resolution and n s is an integer 
that specifies the number of the time step and such that < n s < N r . Each 
Tf- S corresponds to a unique Pk s on the ray path starting from the transmitter 
L such that Pk s can be reached from L in time Tk s - Next, for each time step 
corresponding to time Tk s to reach candidate Pks from L, we step through the 
range of initial angles for the signal starting from the receiver S and check 
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whether the curve from S will intersect Pks in time tk — t& s . Alternatively, we 
can step through a range of reflection angles at the point Pks that is found 
for the given Tk s and check whether the reflected signal will intersect S in time 

tk — T ks ■ 

6 Reconstruction architectures 

When a wave signal reflects from an obstacle it forms a broken ray. Consider 
the broken ray LP {J PS with transmitter at point L, receiver at point S and 
reflection point P. Consider the smallest time period T for which the transmitter 
at L sends beams along all angles from its range of angles. In K 2 this range 
is < 9x,...,6n < 27r. In a high frequency system, T is sufficiently small that 
the shape and position of the obstacles can be considered constant for the time 
period T. The coordinates of each ray can be considered a vector that contains 
the coordinates of its transmitter and receiver, the angle at which the ray was 
sent from the transmitter and the travel time of the ray from the moment it 
was sent from the transmitter to the moment it was received at the receiver. 
For example, 

B = (xi,yi,x r ,y r ,6,t) 

are the coordinates of a ray in M 2 with transmitter (xi,yi), receiver (x r ,y r ), 
angle of transmission 9 and travel time t. In R 3 , 

B = (xi,yi,zi,x r ,y r ,z r ,tf>,0,t) 

are the coordinates of a ray with zenith and azimuth angles of transmission <j) 
and 9. It is important to emphasize that the position of the receiver is not 
known in advance. Only when the ray is received at a receiver the coordinates 
of the receiver are known or measured, and only then the ray coordinates i.e. 
the data point composed of measurements is formed. 

It is possible to design a communication protocol, e.g. a packet based net- 
working protocol, that sends ray coordinates or portions of ray coordinates 
between transmitters and receivers. Modern computer networking is based on 
theory and software for packet based networking protocols and UCLA has pio- 
neered theory and software for packet based communication protocols. 

My design goal for the communications protocol between transmitters and 
receivers is to make it as power efficient, fast and reliable as possible. In order to 
achieve this goal, the protocol and architecture of the communications network 
between transmitters and receivers should enable as much information as pos- 
sible to be inferred rather than communicated. Only the minimum amount of 
information or energy should be sent between transmitters and receivers. The 
network architecture for broken ray communications is composed of transmit- 
ters, receivers and an analyzer or analyzers that combine information from the 
transmitters and receivers and can infer additional information that is not di- 
rectly communicated. The routers in this active networking architecture are the 
obstacles. By virtue of braking the communications ray sent from a transmitter 
an obstacle routes information about its position to the receiver. 

In addition, the receiver or analyzer must infer the time period for which 
a ray was sent. There are several technical approaches that can be used in 
order to achieve this. One approach is to transmit at different frequencies in 
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a block of adjacent time periods. Another approach is to simply leave a time 
gap between time periods so that all sent signals for the preceding period will 
arrive before the signals from the next time period. The receiver or analyzer 
can also infer the angle of transmission. In order to enable this, a transmitter 
can send different angles or combinations of angles at different frequencies and 
in this way each communications frequency implies an angle or combination 
of angles of transmission. The set of unique frequencies £i,£2>— iO use d for 
communication in the network is so important that is should be considered a 
part of the communications architecture and will be referred to as a pallette. In 
addition, the frequency at which a given ray was communicated should be part 
of the ray coordinates: 

B = (xi,yi,zi,x r ,y r ,z r ,<j),6,t,^) 

are the coordinates of a ray communicated at frequency £. 

For each time period T n , the analyzer should then be able to infer a list of 
ray coordinates 

B n ,N = ( x l n , N i yi n , N i z l n ,N J x r n ,N -i Ur n ,N -i z r n ,N ! 071, iV j @n,N j t-n,N i £,n,N ) 

For example, an anazyler that is connected to all transmitters and receivers 
records that a signal is sent with initial angles (f> and 9 from transmitter with 
coordinates xi,yi,zi at time t star t at frequency £. It then records that re- 
ceiver with coordinates x r ,y r , z r has received the signal with frequency £ at 
time t en d- The analyzer computes t = t en d — t start and forms the data point 
B = (xi,y h zi,x r ,y r ,z r ,(f),6,t,£). 

Then by knowing the speed of sound for each frequency, and for ultrasound 
it is the same for different frequencies, and by direct application of formulas 16 
or for variable speed of sound the algorithm from section 7, the analyzer can 
infer the coordinates 

PnS = (ln,l,!/n,l,Zn,l) 
Pn,2 = (x„,2,y n ,2,Z n ,2) 

(48J 

P(n,N) = (Xn,N,y n ,N,Z n ,N) 

of the broken ray reflection points for the given time interval T„. This list of 
coordinates can also be written as 

P(Vi, T n) = (x m ,T n ,y m ,T n ,z„ u T n ) 

(49) 

P(VN,T n ) = (x m ,T„,y VN ,T n ,Z m ,T„) 

and gives a discretization of dfli(t). By the definition of the shape and 
trajectory reconstruction problem, for each time interval T, we are given at 
least one broken ray or data point B that reflects at any given point P on the 
obstacle's boundary. Therefore, we can then reconstruct the whole boundary. 
Section 8 proves that this reconstructed boundary is unique. In a version of the 
shape and trajectory reconstruction problem when fewer data points are given, 
portions of the boundary can still be reconstructed. 
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7 Reconstruction algorithms 



The input to the following algorithm is the speed of sound c(x) for the domain 
f2 and a set of data points or ray coordinates corresponding to broken rays. The 
output is a set of points in R 3 . The points from the output arc the reflection 
points reconstructed from the input data. 

Require: Set of broken ray data points B k = (xi,yi, zi,x r , y r , z r ,<j>k,9k,tk,€k) 

Require: Speed of sound c(x) for domain fi 

{Algorithm for Shape and Trajectory Reconstruction of Moving Obstacles} 
{Estimated time complexity is 0(T 2 A) where T is the number of discretiza- 
tion points for the time of flight, and A is the number of discretization points 
for the angle space} 
for all data points B k do 
hu = 

L = (Xo,Yq, Zq) = (xi,yi,zi) set this initial position to be position of 
transmitter 

S = (aXo,aYo,aZQ) = (x r , y ri z r ) set this initial position to be position of 
receiver 

$o = 4>k 
e = e k 
To = o 

aT = 

for s = N r - 1 do 

{Compute the next point on the ray from the transmitter by fourth order 
Runge-Kutta step and the ray tracing system 42} 
X s+1 = RKA X {h k ,T s ,X s ,Y s ,Z s ,<f> s ,Q s ) 
Y s+1 = RK4 Y (h k , T S1 X S , Y s , Z S ,$ S ,6 S ) 
Z s+1 = RKA z {h k , T s , X s , Y s , Z s , $ s , 6 S ) 
= RKA^{h kl T s , X S ,Y S , Z..<K. (-).:■ 
Q s+1 = RK4 e (h k ,T s ,X s ,Y s ,Z s ,$ s ,O s ) 
T s +i =T S + h k 

P s +i = (X s+ i,Y s+ i, Z s+ i) point on solution of ray tracing equations 
with initial values for transmitter that is at time T s+ \ away from the 
transmitter L 
if !(P s+ i € fi ) then 

There must be a measurment error. Continue with next data point B k 
end if 

for all angles a$o>a©o in discretized angle space do 
for p = ->• N r - 1 do 

{Compute the next point on the ray from the receiver by fourth 

order Runge-Kutta step and the ray tracing system 42} 

aX p+ i = RK4 x {h k ,aT p , aX p , aY p ,aZ p , a$ p , a@ p ) 

aY p+ i = RKAy{h k , aT p , aX p , aY p , aZ p , a$ p , aO p ) 

aZ p+1 = RK4 z (h k ,aT p ,aX p ,aY pi aZ p ,a$ p ,aQ p ) 

a$ p +i = RK4<s,(h k , aT p , aX p , aY p , aZ p , a$ p , a<d p ) 

aO p +i = RK4Q(h k , aT p , aX p , aY p , aZ p , aQ p , a6 p ) 

aT p+ i = aTp + h k 

Pa p+1 = (a>Xp+i, aYp+i, aZ p+ i) point on solution of ray tracing 
equations with initial angles a$o an d a6o and initial position S, 
that is time aT p+ i away from S 
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if !(P Qp+1 G flo) then 

Exit this for loop and continue with next pair of initial angles 

a^OjaQo from outer for loop 
end if 

if distance(P s+ i, P ap+1 ) < ei and \T S+1 + aT p+ i - t k \ < e 2 then 
Pk = P s +i {Solution for current data point Bk found. Continue 
with next data point Bk+i} 

end if 

if T s+ i + aT p+ i > t k + e 2 then 

{We are over the travel time budget tk- Continue looking for a 
solution with the next set of initial angles a$o, a©o-} 
end if 
end for 
end for 
end for 
end for 

RK4 stands for a fourth order Runge-Kutta method although other time- 
dependent numerical methods can be used as well. The above algorithm leads 
to a class of new algorithms when data structures are used to store the rays 
for faster processing. For example, when rays from the receiver in its set of 
directions are stored and looked up in a data structure such as an array, a hash 
table or a set of hash tables, the performance of the algorithm can be improved to 
0(TA) where T is the number of discretization points for the time of flight, and 
A is the number of discretization points for the angle space. These algorithms 
are very suitable for parallelization and their parallelized implementation is key 
for efficient real-time processing and performance. 

When the above algorithm is run on a set of points {Bk} from one sampling 
time interval Tk the algorithm reconstructs the shape of the obstacle during 
this sampling interval. In order to reconstruct the trajectory of the obstacle, 
the algorithm is run on the data points for each of the sampling intervals. 

8 Existence and uniqueness of a reconstructed 
point 

We consider an environment without caustics therefore in order to address the 
question of existence and uniqueness of a reconstructed point we should con- 
sider the case when the receiver is located on the path of an unbroken ray. In 
that case, the above algorithm could return any point of the unbroken ray to 
be a reflection point from an obstacle. Therefore, we should check that the rays 
from the data points that are given as input to the above algorithm do not 
correspond to unbroken rays. One straightforward approach is to trace the ray 
from the transmitter for the whole travel time and then check whether it inter- 
sects the receiver from the same data point. Broken rays are not diffcrcntiablc 
at the reflection point and unbroken rays are differentiable, therefore, if c{x) 
is differentiable and the traced ray from the transmitter intersects the receiver, 
Fermat's principle of least time implies that there is no broken ray starting from 
the same transmitter and with a reflection point on the unbroken ray that has 
the same travel time and that is received by the same receiver. Alternatively, 
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we can precompute and store in a table the initial angles and travel times for 
unbroken rays starting from a transmitter and ending at a receiver and filter 
out rays with such initial angles from the input to the algorithm from section 7. 
The cost of this filtering step is O(T) for each data point where T is the number 
of discretization points for the travel time. The implementation of this filtering 
step can be incorporated in the above algorithm: if the ray from the transmitter 
is traced to the end of the time budget and the end of this ray intersects the 
receiver then the current data point is implicitly filtered out because no unique 
solution for the reflection point is reconstructed. If the ray from the transmitter 
is traced for the travel time, and this ray does not intersect the receiver, and 
no reflection point is found, then its data point is filtered out because it is a 
measurement error and no solution for this data point exists. Therefore, when 
we have a set of data points for broken rays without measurment errors such 
that this set satisfies the conditions of the shape and trajectory reconstruction 
problem, a unique solution to this problem exists and can be found by providing 
the set of broken rays as input to the algorithm from section 7. 

9 Experimental results 

This section presents results from numerical software experiments. We consider 
a basic example in which a small circular reflecting obstacle is moving in the xy 
plane on a circle of radius 1 with center at the origin of our Cartesian coordinate 
system. Let the speed of sound be 1. Position a transmitter and a receiver at 
the origin. Then the travel time of a reflected signal will be 2. Consider the 
moment when the obstacle is at point {l/^/2 1 l/\/2). When the transmitter is 
sending signals with initial angles 8 such that 8 ^ j i.e. \0 — j \ > eg then no 
reflected signal will be detected by the receiver. When the transmitter sends a 
ray at initial angle of 8 = j, the reflected signal will be received by the receiver. 
This gives the data point B k = (0, 0, 0, 0, 0, 0, § , f , 2, £ fc ). We deduce a = 1, 
c = and by 16 reconstruct the reflection point 

x k = l/V2,y k = l/V2, z k = 

In order to reconstruct the shape of the obstacle, we would like to reconstruct 
a neighborhood of a point. This is accomplished by a distributed system of 
transmitters and receivers and by varying the transmission angles. We place a 
second transmitter and a second receiver at the point (-\/2, 0, 0). By sending a 
ray at an initial angle 8 close to j we reconstruct a point close to (x k = ^iVk = 
■^, z k = 0). We reconstruct more points on the surface of the obstacle by 
varying, for example by rotating the beams, the initial angles of the transmitted 
rays as well as by varying the positions of the transmitters and the receivers, 
for example by placing transmitters and receivers at different positions around 
the obstacle. By placing transmitters and receivers in certain configurations, 
for example only in front of the obstacle, we can recover portions of the shape 
of the obstacle. In order to reconstruct the whole boundary, we have to choose 
transmitter and receiver configurations that allow the collection of data for 
each point on the obstacle's boundary. We repeat the reconstruction procedure 
for different sampling intervals during which the obstacle can be considered 
stationary and recover the trajectory of the obstacle. A sampling interval is a 
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short time period during which the obstacle is considered stationary and during 
which we collect data points Bk for reconstructing the shape of the obstacle 
during this time frame by varying the initial angles of transmission as well as 
the positions of the transmitters and the receivers. For example, at the next 
sampling interval, the obstacle can be at position (0, 1). When the transmitter 
at the origin sends a ray at angles = | , <j) = § the receiver will receive a signal 
with travel time 2. We deduce, a = 1, c = and by 16 the reflection point (0,1). 

I have developed a multi-threaded Java implementation of the reconstruction 
algorithm for variable speed of sound from section 7 and numerical results are 
presented in the following tables. Table 1 shows reconstructed trajectory of a 
point moving on the part of the unit circle that is in the first quadrant. The 
first eight columns of the table contain the reconstruction data or input for the 
algorithm and the last three columns contain the x, y and z coordinates of the 
reconstructed points. 
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zl 


xr 


yr 
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xp 


yp 


zp 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.00 


0.98 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.06 


0.98 


0.06 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.13 


0.99 


0.13 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.19 


0.96 


0.18 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.25 


0.97 


0.25 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.31 


0.93 


0.30 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.38 


0.93 


0.37 


0.00 



Table 1: Reconstruction of a point moving with constant speed c=l in the first 
quadrant on an arc of the unit circle for a fixed signal frequency £ and constant 
travel time T=2. 

Consider a circular reflecting obstacle in the plane xy moving away from the 
origin along the line x=y in a medium with variable speed of sound c(x, y) = 
x + y + 1. We place a transmitter and a receiver at the origin. In this case, the 
domain Qq is a circle of sufficiently large radius that contains the origin. Table 
2 shows the trajectory of a point on the obstacle on the line x=y corresponding 
to data with different travel times. 
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xr 


yr 
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T 


xp 


yp 


zp 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


2 


1.55 


1.55 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


3 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


4 


7.89 


7.89 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


5 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


6 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


7 


68.37 


68.37 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.57 


0.79 


8 


138.15 


138.15 


0.00 



Table 2: Reconstruction of a point moving with speed c(x, y) = x + y + 1 on 
the line x=y for a fixed signal frequency £. The initial transmission angles are 
4> = f and 6 = f 

We check whether the computation of the above table by the Java program 
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is correct as follows. The time for the ray to reach to obstacle can be computed 
by the formula 



Therefore, 




By symmetry, for this particular example, this time t is half of the total 
travel time T. Then for a travel time T = 2, or t = 1, we compute X = Y = 
1.55. This result matches the corresponding result for px and py from Table 
2 obtained by numerical integration. In order to reconstruct a neighborhood 
of a point, we place more receivers and vary the transmission angles or the 
location of the transmitter. For example, we place adjacent receivers on the x 
axis covering the interval [-N,N], where N is a constant that depends on the size 
of the neighborhood that we would like to reconstruct, and also send signals 
with initial = ? ± e. 

10 Design considerations 

One classical approach to computing the time of flight is to use synchronized 
clocks. There are elegant protocols for network time synchronization and a lot 
of computing and communications devices already have synchronized clocks. In 
other scenarios particularly in environments with resource constraints and high 
reliability requirements, the need to maintain a clock creates an additional point 
of failure and adds cost and complexity to the system. 

Another well-known approach to computing the time flight is to use the 
phase shift of the signal and potentially the use of more than one frequency in 
order to compute the time of flight reliably This approach is more useful for 
high speed signals. 

One can consider a model of the transmitter to be a system that can rapidly 
send rays at dynamically programmable zenith and azimuth angles. One exam- 
ple of such a system are ultrasonic phased arrays. Ultrasonic phased arrays are 
a well-known and mature technology and offer an off-the-shelf component for 
the implementation of the transmitter device for tomography in the presence of 
obstacles [L] as well as for shape and trajectory tracking of moving obstacles. 

One can consider a model of the receiver to be a sensor array that can convert 
incoming wave signals into signals that can be analyzed by a computer system. 
In the case of ultrasonic signals, a receiver is an ultrasound sensor. Ultrasound 
sensors are a well-known and mature technology that is reliable, efficient and 
cost-effective. 

11 Reconstruction algorithms for infinite travel 
times and absorbing obstacles 

Consider a set of synchronized transmitters that send signals or rays at a range 
of transmission angles. We consider the case when these rays are not received 
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by any receiver and define each such ray by the coordinates of its transmitter 
and the angle of transmission. We refer to such rays as lost rays and consider 
the physical model when a lost ray is absorbed by an obstacle. The coordinates 
of a lost ray B in R 2 are defined as 

B = (x h y h 0,O 

and in R 3 

B = (xi,yi,zi,cf),6,€) 

where <j> and are the zenith and azimuth angles of transmission, and xi, yi, 
Zi are the coordinates of the transmitter endpoint of the ray and £ the signal 
frequency. 

The version of the shape and trajectory reconstruction problem for absorbing 
obstacles and infinite travel times can be defined in terms of sets of lost rays: 
Consider a bounded convex domain fio € R 3 with one moving obstacle fti(t) C 
fioi where tel. Let = fi \fii(*) where fii(t) is a set in M 3 . Let C 
fio. Suppose that for all t we are given a set T of lost rays in fl(t) such that 
for each point P <G dQ.\{t) there is at least one lost ray 7 <G T absorbed at P. 
The shape and trajectory reconstruction problem for an absorbing obstacle is 
to reconstruct dfli(t) given the sets T(t). 

Not all rays will be lost and the receivers will receive signals that are not 
absorbed by absorbing obstacles. The coordinates of these unbroken rays are 
defined to include the coordinates of the receiver. In R 2 the ray coordinates are 
defined as 

B = (xi,yi,x r ,y r ,6,£) 

and in R 3 

B = (xi,yi,Zl,X r ,y r ,Zr,4>,Q,Q 

where 4> and are the zenith and azimuth angles of transmission, xi, yi, zi are 
the coordinates of the transmitter endpoint of the ray, and x r , y r , z r are the 
coordinates of the receiver endpoint of the ray and £ is the signal frequency. 

Consider the smallest time period T for which a transmitter at point L 
sends beams along all angles from its range of angles. In R 2 this range is 
< 0i, On < 27r. In a high frequency system, T is sufficiently small so that 
the shape and position of the obstacles can be considered constant for the time 
period T. 

For each time period T, an analyzer such as a computer system connected 
to the transmitters and receivers, gets as its inputs the list of all coordinates 
of all lost and unbroken rays. Then for the given inputs, the analyzer runs the 
following algorithm for reconstructing dfl\ for the given time period T: 
Require: Set of all image points in Hq and a set of lost and unbroken rays B 
Require: Speed of sound c(x) in fi 

{Estimated time complexity is O (number of image points in fio)} 

Paint all image points p € fio gray 

for all ray B trace ray B given its initial conditions until the ray reaches dfta 
do 

for all point p € B such that p € £Iq do 
if B is unbroken then 

Paint p white 
else if p is gray then 
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Paint p black 
end if 
end for 
end for 

The output of the above algorithm for absorbing obstacles will be an image 
of f2o that paints the boundary of the absorbing obstacles in black, the medium 
n in white, and points from Qq that were not touched by any of the input rays 
in gray. 

The initial value ray tracing in the above algorithm is done by the system 42 
and fourth order Rungc-Kutta time steps as in the algorithm from section 7. In 
the above algorithm for absorbing obstacles the ray tracing continues until the 
ray reaches dQo while in the algorithm for reflecting obstacles from section 7 
additional conditions are that the ray tracing continues only for the duration of 
the time budget or until a reflection point is found. The computational efficiency 
of the algorithm for reflecting obstacles is significantly better because it does 
not require processing on the order of all points in the domain. It is faster to 
track reflecting obstacles with the algorithms from this work. 

In the above algorithm for absorbing obstacles, we do not make any assump- 
tions on how the rays are generated and scan JIq. We will consider next two 
scanning algorithms that will in turn lead to new imaging algorithms. First, 
consider a three dimensional scanning algorithm in which the ray sweeps across 
a plane. Then if for all time intervals T of the plane scan the receivers receive 
a signal from the transmitter, the plane or angle sector is marked as a white 
plane. If for all time intervals T of the plane scan the receivers do not receive 
a signal from the transmitter, the plane or angle sector is marked as a black 
plane. When an individual signal is missed, the data is recorded as a black ray. 
When an individual signal is recorded it is marked as a white ray. Therefore, 
the input of the new imaging algorithm is a series of black and white rays and 
planes or plane regions between pairs of rays, fio is then painted using planes 
or sectors in addition to rays. 

11.1 Reconstruction algorithms for absorbing and reflect- 
ing obstacles 

We combine the approaches from this section and section 7 into a two-phase 
algorithm to reconstruct the boundary and trajectory of one reflecting obstacle 
and many absorbing obstacles or one obstacle that has both absorbing and 
reflecting regions on its boundary. For a given sampling period T: 
Require: Set of image points in fio and a set of rays B 
Require: Speed of sound c(x) in fio 

{Estimated time complexity is O (number of image points in fio)} 

{Phase 1} 

{Reconstruct boundary of reflecting obstacles} 
Paint all image points p £ flo gray 
for all broken rays B do 

Reconstruct reflection point of B by reconstruction algorithm for finite 

travel times and reflecting objects 

Paint this reflection boundary point red 
end for 
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{Phase 2} 

{Reconstruct boundary of absorbing obstacles} 
for all ray B do 

for all point p € B such that p € Hq do 

if B is broken or unbroken and p is not black or red then 

Paint p white 
else if p is gray then 

Paint p black 
end if 
end for 
end for 

The output of the above algorithm for an absorbing and reflecting obstacle 
will be an image of fio that paints the boundary of the reflecting obstacle in 
red, the boundary of the absorbing obstacle in black, the medium £1 in white, 
and points from Qq that were not touched by any of the input rays in gray. 

12 Applications 

The theory for shape and trajectory tracking of moving obstacles has broad 
applicability. In this section, I present an application of shape and trajectory 
tracking of moving obstacles to Internet telephony and mobile communications. 
Shape and trajectory tracking of moving obstacles that reflect ray signals sent 
from moving Internet telephone transmitters and receivers enables efficient one- 
hop routing of messages through reflection from mobile telephone transmitters 
and receivers to the nearest Internet router or peer. 

A router in this architecture is a transmitter and receiver of ray signals. 
Moreover, every every mobile phone communicator can participate in a peer-to- 
peer network where every transmitter and receiver is a router and/or reflector. 
This mobile Internet phone architecture is illustrated in Fig. 2. 

The combination of the methods for obstacle tracking from this work with 
the numerical methods for tomography in the presence of obstacles from [L], 
leads to a new method for direct mobile communications with reflected rays 
based on tracking the location of reflectors and finding the speed of sound in 
the environment. Through reflection, messages can be communicated around 
obstacles blocking the line of sight as well as around areas with a high error 
rate. The reflectors for the shape and trajectory tracking method do not have 
to be installed in order to enable a network and can be obstacles such as people, 
buildings, trees or cars. The two initial angles provide a way of addressing the 
obstacle and the receiver i.e. a way of calling the obstacle and the receiver by 
dialing two angles of transmission. When there are caustics or blocks on the 
path from the transmitter K to the receiver R then K and R can communicate 
through an obstacle/reflector RK. 

12.1 Optimal communication path 

In order to increase reliability or find a faster/better route, we can send the 
message along several different angles to several different reflectors. 

We define a new physical layer protocol for point to point communication 
called IHOP. The goal of this protocol is to route messages in one hop. The 
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Mobile 
Internet 
Phone 4 




Figure 2: Mobile Internet Phone Architecture. Mobile Internet Phone 1 trans- 
mits a message to Mobile Internet Phone 2 via the broken ray Ray 1. Knowing 
the speed of sound in the environment, Mobile Internet Phone 2 applies the 
algorithm for shape and trajectory tracking of moving obstacles to reconstruct 
the reflection point at Obstacle 1 and send back response signals in the direc- 
tion of this reflection point. Note that the algorithm from section 7 enables the 
reconstruction not only of the reflection point but also of the arrival angles of 
the signal at Mobile Internet Phone 2. The fastest/most efficient reconstruction 
is when communication is at relatively short distances and the speed of sound 
is approximately constant and the ray paths are approximately linear segments 
as illustrated by Ray 2 and Ray 3. 
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IHOP address of a target host is always relative to the position of another host 
called the origin and consists of this position and a pair of zenith and azimuth 
angles or direction cosines for reaching the target host by sending a signal from 
the origin to the target. The signal may reach the target by an unbroken ray 
or by a broken ray when it reflects from an obstacle. 

For improved bandwidth and reliability, we can reach the host over different 
angles. The IHOP address of the host specifies a ray. 
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